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DOCUMENTATION  FOR  PCDIST,  A  FORTRAN  SUBROUTINE 
TO  CALCULATE  POLLACK-CUZZI  SCATTERING 


1 .  INTRODUCTION 

Hodkinson  and  Greenleaves'  describe  a  method,  based  on 
classical  diffraction  and  geometric  optics,  for  calculating  light 
scattering  from  a  polydispersion  of  transparent  spheres.  The 
scattered  field  is  treated  as  arising  from  three  independent 
processes:  diffraction  of  light  rays  passing  near  the  spheres; 

reflection  of  rays  incident  on  the  spheres;  and  refraction  of 
rays  passing  through  the  spheres.  The  three  components  are 
summed  without  regard  to  phase  relations  on  the  assumption  that 
the  size-range  of  spheres  is  sufficiently  wide  to  smear  out  the 
far  field  interference  effects  created  within  any  particular 
sphere.  An  exact  calculation  of  the  scattering  by  a  polydis¬ 
persion  of  spheres  can  be  carried  out  using  Mie  theory  (in 
comparisons,  the  approximation  of  Hodkinson  and  Greenleaves  holds 
up  well) ,  but  the  computation,  while  commonplace  today,  was 
challenging  for  computers  of  the  time,  and  few  tables  of  Mie 
results  were  available  then.  The  purpose  of  Hodkinson  and 
Greenleaves'  work  was  to  provide  a  method  for  solving  many 
practical  scattering  problems  while  avoiding  the  use  of  Mie 
theory . 


Contemporary  efforts  to  calculate  light  scattering 
from  irregular  particles  has  much  in  common  with  the  earlier 
sphere  scattering  problem.  The  detailed  structure  of  the 
scattering  pattern  produced  by  any  given  arbitrary  particle  in  a 

particular  orientation  is  very  difficult  today  -  usually 

impossible  -  to  calculate;  in  the  ordinary  circumstance  of 

scattering  by  a  cloud  of  many  rough  particles  in  random  orien¬ 
tations,  the  fine  details  of  the  picture  should  all  be  averaged 
out.  So  rather  than  attempting  a  (probably  hopeless)  detailed 
calculation  for  cvnry  irregular  particle  in  an  ensemble  and  then 
averaging  the  results,  one  would  like  to  develop  and  apply  a 
relatively  simple  procedure  that  can  give  the  average  scattering 
pattern  directly. 

Pollack  and  Cuzzi^  suggest  using  the  same  three- 
component  approach  to  nonsphere  scattering  as  had  been  success¬ 
fully  employed  by  Hodkinson  and  Greenleaves  to  simplify  scat¬ 
tering  from  spheres.  They  parameterized  particle  shape  in  terms 
of  sphericity  and  surface  roughness,  and  wrote  down  rules  for 
calculating  and  combining  the  diffraction,  reflection,  and 
refraction  components  of  scattering.  This  is  an  empirical 
approach,  and  it  succeeds  according  to  whether  calculated  results 
and  light  scattering  measurements  agree  with  each  other.  Pollack 
and  Cuzzi  examined  experimental  light  scattering  measurements, 
previously  published  by  a  number  of  authors,  and  found  that  their 
formulae,  with  appropriately  chosen  parameter  values,  could  fit 
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the  experimentally  determined  phase  functions  of  irregular 
particle  systems  very  well. 

To  make  the  Pollack  and  Cuzzi  (P-C)  approximation 
readily  available,  I  have  written  a  Fortran  subroutine  that  can 
be  called  by  a  main  program  to  which  the  subroutine  returns  the 
P-C  phase  function  and  average  cross-sections  for  a  distribution 
of  irregular  particles  specified  by  a  "parameter  array"  in  the 
main  program.  This  report  is  the  documentation  for  that  sub¬ 
routine,  PCDIST.  I  shall  review  the  P-C  approximation  more 
fully,  describe  how  it  is  implemented  in  PCDIST,  and  give  a  few 
examples  of  its  use. 


2.  POLLACK  AND  CUZZI 'S  SEMI-EMPIRICAL  THEORY 

Pollack  and  Cuzzi 's  theory  treats  the  scalar  (unpolar¬ 
ized)  scattering  properties  of  groups  of  randomly  oriented 
nonspherical  particles.  They  begin  with  a  size  parameter 
distribution  n(X)  of  particles,  where  the  size  parameter  of  a 
nonspherical  particle  is  equal  to  the  size  parameter  (circum¬ 
ference  to  wavelength  ratio)  of  a  sphere  with  equal  volume.  That 
part  of  the  distribution  with  size  parameters  less  than  or  equal 
to  some  breakpoint  XO  is  treated  as  if  it  consisted  of  perfect 
spheres,  and  Mie  theory  is  applied  to  determine  the  phase 
function  and  the  efficiencies  for  absorption,  scattering,  and 
extinction.  Larger  particles  are  handled  using  the  three- 
component  approximation. 

The  diffraction  component  of  the  scattered  light 
produced  by  an  irregular  particle  of  size  parameter  X  >  XO  is 
taken  to  be  the  same  as  that  produced  by  a  sphere  of  size 
parameter  /RX,  where  the  factor  /R  (R  is  the  ratio  of  the 
particle's  surface  area  to  the  equi-volume  sphere's  surface  area) 
accounts  for  the  irregular  particle's  greater  (on  average) 
projected  area.  In  general,  the  diffracted  intensity  distribu¬ 
tion  about  a  sphere  of  size  parameter  x  is  given  by  equation  1.’ 


Id{9,x)  =  Cd 


2Ji(x  sin0) 
X  sind 


(1) 


where 


k—  ^{1  +  cos^  B). 


(2) 
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Cp  is  a  normalization  constant,  chosen  so  that 


dSl  =  2Tr 


Id 


sin  9  d0  =  4ir. 


(3) 


Pollack  and  Cuzzi  assume  that  the  geometrical  optics 
expression  for  external  reflection  from  large  spheres  applies 
c‘iSO  to  all  irregular  particles  with  X  >  XO.  That  expression  is 


r  siD(f)  -  [|m|^  -  1  -hsin^(f)l^ 

^  2  ^\sm(|)  +  [|mP  -  1  +  sin^(®)]U 

-l-ic  f 

2  |m|2sin(f)  +  (|mp  -  1  +sin^(|))5  / 


where  m  is  the  complex  refractive  index,  and  is  a  normal¬ 
ization  constant  so  that  Ig  also  satisfies  equation  3.  The  two 
terms  represent  reflection  of  parallel  and  perpendicular  polar¬ 
izations,  respectively. 

These  first  two  components  are  the  same  as  those 
Hodkinson  and  Greenleaves  used  in  their  approximation  for  scat¬ 
tering  from  dielectric  spheres.  Irregular  particle  scattering  is 
distinguished  from  sphere  scattering  in  Pollack  and  Cuzzi 's 
treatment  of  the  transmitted  (i.e.,  refracted)  component  of  the 
scattered  light.  In  their  scheme,  the  angular  redistribution  of 
light  having  entered  the  particle  is  governed  by  an  empirical 
law. 


If  =  07'exp(l  +  bO) 


(5) 


where,  again,  Cj-  is  a  normalization  constant.  The  empirical 
parameter  b  can  be  related  to  another  constant,  G,  by 


_  fo^  lr  M 
-  i’.hdB’ 


(6) 


which  relates  the  energy  reradiated  into  the  forward  hemisphere 
to  that  reradiated  into  the  back  hemisphere.  For  smooth  parti¬ 
cles  that  do  not  deviate  the  incident  light  very  much,  G  is 
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large,  -10;  whereas,  for  very  rough  particles  that  redistribute 
the  light  more  uniformly,  G  is  small,  -1.  Pollack  and  Cuzzi  do 
not  suggest  any  quantitative  means  for  assigning  a  value  to  G 
from  an  examination  of  the  irregular  particles;  in  fact,  G  is 
principally  used  as  a  fudge-factor  to  bring  about  agreement 
between  measured  and  calculated  phase  functions. 

The  scattering  efficiency  associated  with  each  of  the 
three  components  must  be  specified  so  they  may  be  properly 
weighted  when  summed.  Consider  a  single  irregular  particle  (or 
infinitesimally  narrow  range)  of  size  parameter  X  (>X0) .  The 
diffraction  efficiency  is  taken  to  be  exactly  unity  (Q^  =  1). 

That  is,  it  is  assumed  that  any  large  irregular  particle 
diffracts  an  amount  of  energy  equal  to  that  falling  on  its 
(average)  physical  cross-section.  The  efficiency  factor  for 
external  reflection  is  =  {l/C^}  where  is  the  normalization 
factor  in  equation  4.  Finally,  the  scattering  efficiency  for  the 
transmitted  component  is  given  by 


Qt  =  Qs  -  Qd  -  Qn 


(7) 


where  is  the  Mie  scattering  efficiency  for  the  equivalent 
volume  sphere.  The  total  scattering  efficiency  for  the  large 
particle  is  taken  to  be  RQ^  (because  the  total  scattering  cross- 
section  should  be  proportional  to  the  irregular  particle's 
projected  surface  area) .  The  absorption  efficiency  for  a  large 
particle  is  just  Q^,  the  Mie  absorption  efficiency  for  an  equi- 
volume  sphere,  unless  the  particle  is  highly  absorbing  (2jnX  >  1)  , 
in  which  case,  an  extra  factor  R  is  also  included. 

The  composite  normalized  phase  function  for  a  large 
irregular  particle  is  given  by 


p  r  _  r  Q  P  ,  r  Q  R  ,  T  Qt 


(8) 


3.  SUBROUTINE  POSING  AND  ITS  DAUGHTERS 

The  Fortran  program  PCDIST  that  implements  the  P-C 
approximation  is  diagramed  in  Figure  1.  The  name  PCDIST  is  used 
to  refer  both  to  the  entire  collection  of  subroutines  shown  in 
the  figure  and  to  the  uppermost  level  subroutine  that  interacts 
directly  with  the  user's  main  program.  The  program  structure  is 
best  explained  by  starting  in  the  middle,  with  POSING. 
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USER’S  MAIN 
PROGRAM 


Establish 
parameter  array. 


Use  or  display 
results. 


Figure  1.  Relationships  Among  Subroutines  Comprising  PCDIST 


11 


PCSING  controls  the  calculation  of  the  normalized 
phase  function  and  scattering  and  absorption  efficiencies  for  a 
single  particle.  The  required  parameters  (X,  CMR,  CM! ,  R,  XO, 
and  G)  are  passed  in  the  argument  list  of  the  calling  statement 
within  PCDIST.  The  values  CMR  and  CMI  are  the  real  and  imaginary 
parts  of  the  particle  refractive  index;  the  other  parameters  were 
defined  in  the  previous  section.  Also  passed  are  DELTHA,  the 
increment  between  scattering  angles  at  which  the  phase  function 
is  to  be  evaluated,  and  the  closely  related  NANG,  the  number  of 
angles  between  0  and  180°,  inclusive,  at  which  the  phase  function 
is  to  be  evaluated. 

The  particle  size  parameter  X  is  compared  to  the 
critical  size  XO;  if  X  is  smaller  then  the  subroutines,  QMIE  and 
PFMIE  are  called,  in  turn,  to  calculate  the  particle's  efficien¬ 
cies  and  phase  function  as  if  it  were  a  perfect  sphere.  Both 
subroutines  are  adapted  from  the  excellent  programs  published  by 
Barber  and  Hill."*  QMIE  comes  from  their  program  SI,  and  PFMIE  is 
adapted  from  their  program  S3.  Both  of  these  Mi.-  subroutines 
call  a  number  of  other  subroutines  as  they  proceed;  all  are 
included  in  the  PCDIST  package. 

If  X  is  greater  than  XO,  QMIE  is  still  called  to 
obtain  the  Mie  efficiencies;  but,  the  phase  function  is  calcu¬ 
lated  with  the  help  of  three  new  subroutines.  PFDIFF  is  called 
for  the  diffraction  part  of  the  phase  function,  as  in  equation  1. 
The  quantity  [J, (z)/z]  is  approximated,  very  accurately  for  all 
z,  using  the  algorithm  by  Press*  for  J,(z) .  Actually,  the 
diffracted  light  expressed  in  equation  1  is  symmetric  around  90® 
and  peaks  equally  strongly  in  the  backward  and  forward  direc¬ 
tions,  which  is  nonsense.  Therefore,  I  cut  the  diffracted 
component  off  after  90®.  That  cutoff  may  put  a  small  nonphysical 
kink  in  the  phase  function  curve  at  90®,  depending  on  the 
intensity  of  the  diffraction  component  at  cutoff.  To  avoid  that, 
I  have  added  an  extra  cos (6)  factor  to  the  expression  for  It 

has  no  significant  effect  at  small  scattering  angles  but  insures 
that  the  diffracted  component  goes  smoothly  to  0  at  90®.  PFDIFF 
returns  to  PCSING  a  1-d  array  of  numbers  (PFD) ,  which  are  the 
diffraction  phase  function  values  (except  for  the  normalization 
constant)  at  all  desired  angles  from  0  to  90®,  and  it  returns  the 
normalization  constant,  called  CDD  here.  The  CDD  constant  could 
be  multiplied  against  each  number  in  PFD  to  return  the  normalized 
phase  function;  but,  to  save  some  time,  that  multiplication  is 
performed  later  in  PCSING  in  conjunction  with  Q-weighting  and 
summing  the  three  components. 

The  reflection  component  is  calculated  with  PFREFL, 
which  proceeds  exactly  as  specified  by  equation  4.  It  returns 
the  array  PFR,  with  unnormalized  reflection  phase  function 
components  from  0  to  180®  in  steps  of  DELTHA,  and  the  normal¬ 
ization  constant  CDR.  Notice  that  the  returned  values  depend 
only  on  the  particle  refractive  index  and  no  other  parameters. 
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Frequently,  it  should  be  the  case  that  all  of  the  particles  in  a 
distribution  of  interest  have  the  same  refractive  index,  and  a 
lot  of  computer  time  could  be  wasted  recalculating  exactly  the 
same  reflection  result  each  time  a  new  particle  from  the 
distribution  is  handed  to  PCSING.  Therefore,  we  use  the  Fortran 
SAVE  statement  to  cause  PCSING  to  keep  copies  of  PFR  and  CDR  in 
static  storage,  along  with  the  refractive  index  values,  renamed 
CMROLD  and  CMIOLD.  Each  time  PCSING  is  called,  the  new  values 
CMR  and  CMI  are  checked  against  the  old  saved  ones;  if  there  is 
no  difference,  PFREFL  is  not  called.  Instead,  the  previously 
calculated  PFR  and  CDR  are  simply  reused. 

The  transmitted  component  is  evaluated  in  the 
subroutine  PFTRAN.  Because  of  the  simple  form  assumed  for  Ij^, 
the  normalization  factor  in  this  case  can  be  written  analyt¬ 
ically.  The  condition 


f  /p  dil  =  27r  /  Cr  exp(l  +  W)sind  =  4ir  (9) 

4ir  Jo 


yields 


Ct  = 


2 

e  Jq  sin  9  d9 


(10) 


The  integral  in  the  denominator  can  be  looked  up.* 
After  some  algebra,  we  may  write  for  the  normalized  transmitted 
component 


It  = 


2(&2  +  1)  M 

— rz - e 

+  1 


(11) 


However,  it  is  G,  not  b,  which  is  the  user-specified 
parameter.  The  relationship  between  them  is  discovered  by 
substituting  equation  5  into  equation  6.  One  finds  then  that 


-21dG 

TT 


(12) 
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Subroutine  PFTRAN  uses  equations  12  and  11  to  return  a  1-d  array, 
PFT,  containing  the  normalized  phase  function  at  angles  between 
0  and  180“  in  steps  of  DELTHA. 

Like  PFREFL,  PFTRAN  depends  only  on  a  parameter  (G) , 
which  may  be  constant  for  the  entire  distribution  of  particles. 

So  again,  to  avoid  unnecessary  calculations,  a  copy  of  PFT  and 
the  old  value  of  G  is  saved  in  PCSING,  and  PFT  is  recalculated 
for  new  particles  only  if  the  value  of  G  is  changed. 

PCSING  assembles  the  final  normalized  phase  function 
(PFS)  array  by  calculating  for  each  angle,  indexed  by  I, 


PFSiJ)  = 


CDD 

QSCA 


PFD{1)  +  QR 


CDR 

QSCA 


PFR{I)  + 


QT 

QSCA 


PFT(I) 


(13) 


where  QSCA  is  the  Mie  scattering  efficiency  returned  by  QMIE. 
Finally,  the  particle's  total  scattering  efficiency  is  given 
as  /R  X  QSCA,  and  its  absorption  efficiency  is  either  QABS 
(returned  by  QMIE)  or  x  QABS,  depending  on  2X  CMI , 

It  is  easy  to  temporarily  modify  PCDIST  so  that  it 
writes  a  file  to  disk  containing  the  normalized  total  phase 
function  and  the  three  (normalized)  partial  phase  functions  that 
went  into  it  as  functions  of  scattering  angle.  A  plot  of 
those  data  helps  one  to  visualize  the  P-C  approximation.  For 
example,  Figure  2  shows  such  a  plot  for  a  rela-tively  small 
rough  particle.  Its  P-C  parameters  were:  X  =  8,  CMR  =  1.5, 

CMI  =  0.05,  XO  =  3,  R  =  1.4,  and  G  -  2.  Figure  3  shows  a 
particle  of  the  same  shape  and  material  but  larger  with  a 
smoother  surface;  X  =  30  and  G  =  8,  with  the  other  parameters  as 
shown  in  Figure  2.  Notice  the  higher  frequency  oscillations  and 
the  steeper  transmission  slope  in  Figure  3  compared  to  Figure  2 . 


4.  SUBROUTINE  PCDIST  AND  THE  PARAMETER  ARRAY 

PCDIST  is  the  top  level  of  the  subroutine;  it  gets 
information  about  the  distribution  of  irregular  particles  from 
the  user's  main  program,  calls  PCSING  repeatedly  to  get  the 
scattering  properties  of  each  representative  particle  of  the 
distribution,  and  combines  the  accumulated  results  to  construct 
the  scattering  properties  (phase  function  and  cross-sections)  of 
the  distribution. 

Data  describing  the  cloud  of  irregular  particles  is 
transferred  by  way  of  a  parameter  array,  which  the  user  must 
construct  in  his  main  program;  the  parameter  array  name  is  passed 


14 


Scattering  Angle 

2.  P-C  Phase  Function  when  X  =  8  and  G  =  2 


Scattering  Angle 

3.  P-C  Phase  Function  when  X  =  30  and  G  =  8 
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to  PCDIST  in  the  calling  statement  argument  list.  The  array  must 
have  seven  columns  and  one  or  more  rows.  Each  row  either  repre¬ 
sents  one  particle  or  one  type  of  particle  in  the  distribution. 
The  first  six  columns  contain  the  Pollack-Cuzzi  parameters 

already  discussed  -  in  order:  X,  CMR,  CMI ,  XO,  R,  and  G.  The 

seventh  column  value  tells  what  fraction  of  the  cloud's  particles 
are  represented  by  the  particle  described  in  that  row.  Usually, 
the  seventh  column  is  just  the  particle  size  distribution 
function  evaluated  at  the  row's  X  value;  however,  sometimes  more 
complicated  expressions  will  be  appropriate.  The  parameter  array 
is  a  powerful  and  flexible  structure  that  allows  one  to  solve 
almost  any  scattering  problem  of  this  type  with  a  few  simple 
lines  of  code. 

When  PCDIST  is  called,  the  argument  list  must  also 
contain  NROWS,  the  number  of  rows  in  the  parameter  array,  and 
DELTHA,  the  step  size  in  degrees  at  which  the  phase  function  is 
desired.  Typically,  about  100  rows  should  suffice  for  specifying 
a  particle  distribution,  in  the  sense  that  a  higher  resolution 
specification  using  more  rows  will  not  yield  a  significantly 
different  result.  The  current  version  of  PCDIST  is  dimensioned 
for  a  maximum  of  5000  rows,  which  should  cover  any  situation. 

That  dimension  can  be  easily  changed;  if  it  is  made  too  much 
larger,  a  user's  program  linked  with  PCDIST  will  not  fit  on  a  DOS 
machine.  The  smallest  value  for  DELTHA  allowed  in  the  current 
version  is  0.05°,  but  it  is  difficult  to  imagine  the  need  for 
such  a  fine  scale.  To  shorten  execution  time,  one  should  set 
DELTHA  no  smaller  than  is  needed  to  get  phase  function  values  at 
whatever  particular  angles  are  desired.  However,  DELTHA  should 
not  be  too  large  either  (>«10°)  because  DELTHA  sets  the 
quadrature  size  in  evaluating  the  phase  function  normalizations, 
and  errors  can  creep  in  if  the  scale  becomes  too  coarse. 

DELTHA  =  1°  will  work  well  in  most  circumstances;  in  any  event, 
DELTHA  must  always  be  chosen  such  that  180/ DELTHA  is  an  integer. 

PCDIST  calls  PCSING  NROWS  times,  getting  back  NROWS 
individual  phase  functions  and  absorption  and  scattering 
efficiencies.  Each  phase  function  is  weighted  according  to  its 
contribution  to  the  total  scattering;  then,  the  functions  are  all 
added  together.  Finally,  the  sum  of  the  weighted  phase  functions 
is  normalized.  The  weighting  factor  for  each  contribution  is  the 
product  of  the  scattering  cross-section  (which  is  proportional  to 
X^  times  the  scattering  efficiency)  times  the  column  seven  factor 
(the  relative  number  of  particles  of  that  type) .  The  weighted 
average  scattering  and  absorption  cross-sections  are  similarly 
calculated . 
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5. 


EXAMPLES 


A  few  examples  of  main  programs  that  call  PCDIST 
should  clarify  how  the  parameter  array  is  to  be  used.  The  first 
is  quite  simple. 

5 . 1  Example  1. 

Calculate  the  phase  function  for  a  cloud  of  irregular 
particles  with  a  log-normal  size  distribution.  The  material  of 
which  the  particles  are  made  has  a  refractive  index  of  1.45  + 
0.05i,  and  the  other  Pollack-Cuzzi  parameters  (valid  for  all  the 
particles)  are  XO  =  3,  R  =  1.3,  and  G  =  5.  The  size  distribution 
peaks  at  X  =  5,  and  the  sigma  is  0.2. 

Figure  4  displays  a  code  to  accomplish  example  1.  The 
parameter  array  (called  PAR  here)  and  phase  function  array  (PF) 
must  always  be  dimensioned  exactly  as  shown  in  line  six  to  agree 
with  their  dimensions  in  PCDIST. 

A  log-normal  distribution  n(X)  is  given  by 


n{X)  cx 


1 

(SG)  X 


(-{\nX-\n  {XM)f\ 

- 2W - ) 


(14) 


where  XM  is  the  geometric  mean  size,  and  SG  is  the  geometric 
standard  deviation.  We  consider  particles  with  sizes  up  to  three 
geometric  standard  deviations  from  the  mean  and  divide  that 
range  into  200  intervals  (lines  19-22) . 

The  parameter  array  PAR,  which  in  this  case  will  use 
only  201  of  the  5000  rows  dimensioned  to  it,  is  set  up  in  lines 
24-33.  The  first  column  of  PAR  holds  the  201  values  of  X,  and 
the  seventh  holds  the  corresponding  probabilities,  from  equa¬ 
tion  14.  It  is  not  necessary  to  normalize  the  probability 
distribution;  numerical  normalization  of  the  final  phase  function 
will  be  done  by  PCDIST  and  will  be  unaffected  by  any  constant 
factors  in  distributions  such  as  equation  14. 

This  program  writes  the  scattering  and  absorption 
cross-sections  to  the  screen  (since  particle  size  parameter  was 
given,  a  wavelength  must  yet  be  specified  to  convert  the  cross- 
sections  into  actual  areas)  and  writes  files  listing  the 
calculated  phase  function  and  the  distribution  function  used 
(equation  14) .  Figure  5  shows  plots  of  the  program  output  files, 
PHASFUNC.DAT  and  DISTFUNC.DAT. 
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PROGRAM 

VERSION 

EX1.FOR 

Mon  Mar  01  08:16:58  1993 

■.|geiv;v 

^  ' 

' ,  S'  ' 

LINE  t 

SOURCE  CODE 

PAGE  1 

1 

PROGRAM  EXl 

2 

c  A  ruin  progrui  to  call  pcdist  for  a  log-nonu) 

3 

c  site  distribution  of  particles. 

5 

REAL*8  XLOW.  XHir-H,  DELX.  XM.  SG.  CMR.  CMI.  XO.  R,  G 

6 

REAL*8  PAR(5000.7),  PF(0:3600),  CSCA,  CASS.  OELTHA 

7 

INTEGER  I,  J.  NROWS,  NANS 

9 

CMR  •  1.45 

10 

CMI  =  0.05 

11 

XO  =  3 

12 

R  =  1.3 

13 

G  =  5 

14 

OELTHA  •  1.0 

15 

NANG  =>  181 

16 

XM  -  5 

17 

SG  -  .2 

18 

19 

XHIGH  =  EXP(  L0G{XH)  +  3.0»SG  ) 

20 

XLOW  =  EXP{  LOG(XM)  -  3.0*SG  ) 

21 

DEU  =■  (XHIGH-XL0W)/2OO. 

22 

NROWS  =  201 

23 

24 

DO  30  1=1, NROWS 

25 

X=XL0W+DELX*REAL(I-1) 

26 

PAR(I.l)  =  X 

27 

PAR(I,2)  =  CMR 

28 

PAR(I.3)  =  CMI 

29 

PAR(I,4)  =  XO 

30 

PAR(I.5)  =  R 

31 

PAR(I,6)  =  G 

32 

PAR(I.7)  =  (l./(SG*X))*EXP(-(L0G(X)-L0G(XM))**2/(2. *50**2)) 

33 

30  CONTINUE 

34 

35 

CALL  PCDIST (NROWS.  PAR,  OELTHA,  CSCA,  CABS.  PF) 

36 

37 

WRITE(*.3)  CSCA 

38 

3  FORMATC  CSCA  =  '  .F8.2,2X, ‘times  lambda  squared') 

39 

WRITE(*.4)  A8S(CABS) 

40 

4  FORMATC  CABS  =  '  ,F8.2.2X. 'times  lambda  squared') 

41 

42 

OPEN (UN I T= 1 1 .  FI LE= ' PHASFUNC . OAT ' ) 

43 

DO  40  J=0, (NANG-1) 

44 

40  WRITE(11.21)  DELTHA*REAL(J),  PF(J) 

45 

21  F0RHAT(F8.2,  E12.4) 

46 

CLOSE (11) 

47 

48 

0PEN(UNIT=12.  FILE-'OI STFUNC . DAT ' ) 

49 

DO  50  1=1. NROWS 

50 

50  WRITE(12.22)  PAR(I.l).  PAR(I.7) 

51 

22  F0RMAT(F8.3.  E12.4) 

52 

CL0SE(12) 

53 

STOP 

54 

END 

Figure  4.  Fortran  Code  to  Solve  Example  1 
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Figure  5.  Plot  of  Data  Generated  by  Running  EXl.FOR 


5.2  Example  2 . 

Suppose  a  cloud  consists  of  uniform  NaCl  cubic 
particles/  2.0  fx  on  edge  (as  might  be  generated  by  a  vibrating 
orifice  aerosol  generator  loaded  with  saltwater) .  What  is  the 
phase  function  of  this  aerosol  when  illuminated  by  an  incan¬ 
descent  lamp  with  a  color  temperature  of  6000  K? 

In  this  case,  the  particles  all  have  the  same  physical 
size,  but  a  distribution  in  the  size  parameter  arises  from  the 
broad  spectrum  of  the  incident  wavelengths.  Several  of  the 
parameters  to  be  loaded  into  the  parameter  array  are  functions  of 
wavelength.  A  sphere  with  diameter 


(15) 


has  the  same  volume  as  a  cube  with  side  a.  The  size  parameter  of 
the  equivalent  volume  sphere  is 


""  ~ir~ 


(16) 


and  the  ratio  of  surface  areas  is 


R  = 


1.2407 


(17) 


The  lamp  spectral  distribution  may  be  taken  from 
Planck’s  Law  which,  for  T  =  6000  K  and  X  in  microns,  has  the  form 


/(A)  a 


1 

A5(eT  -  1) 


(18) 


The  relative  number  of  particles  of  each  size  parameter  (i.e., 
column  seven  in  the  parameter  array)  comes  from  this  expression. 
Almost  all  of  the  lamp's  energy  is  radiated  at  wavelengths 
between  0.2  and  2  n,  which  will  be  taken  as  the  limits  of  the 
spectrum.  We  (rather  arbitrarily)  divide  this  range  into  60 
equal  intervals  and  look  at  61  particle  types. 

Another  wrinkle  in  this  example  is  that  the  refractive 
index  is  different  for  each  particle  type;  it  varies  with  wave¬ 
length.  The  real  refractive  index  for  NaCl  over  the  range  of 
interest  may  be  approximated  (with  X  in  microns)  by 


CMR  =  1.52122  + 


01095 

A2 


.00000683 

A< 


(19) 


obtained  by  fitting  the  Cauchy  dispersion  formula  with  data  from 
the  AIP  handbook.’  The  imaginary  part  is  zero. 

Figure  6  is  a  simple  program  that  loads  a  parameter 
array  in  accordance  with  equations  16-19  and  writes  a  file  of  the 
phase  function  values  every  1®.  In  Figure  7,  the  phase  function 
is  plotted. 

5 . 3  Example  3 . 

In  this  last  example,  the  same  conditions  in  example  2 
hold,  except  that  now  we  also  include  a  distribution  in  the 
physical  sizes  of  the  salt  cubes.  Suppose  the  cube  lengths  are 
normally  distributed,  with  a  mean  size  of  2.0  m  and  a  standard 
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Figure  6.  Fortran  Code  to  Solve  Example  2 
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deviation  of  0.3  n,  from  1.0  to  3.0  m-  The  probability  that  a 
random  cube  has  an  edge  of  length  a  is  proportional  to 


P 


(20) 


We  divide  the  particle  size  range  evenly  from  1.0  to 
3.0  M  in  steps  of  0.04  /i  and  must  now  run  over  all  51  particle 
sizes  at  each  of  the  61  wavelengths.  This  is  done  in  the  code  of 
Figure  8,  where  61  x  51  =  3111  rows  of  particle  types  are  defined 
with  a  nested  do-loop.  The  weight  given  to  each  type  (column 
seven)  is  just  the  product  of  the  wavelength  distribution  and  the 
particle  size  distribution  functions.  As  it  happens,  the  phase 
function  calculated  for  example  3  is  only  very  slightly  different 
from  that  of  example  2  because  smoothing  over  size  parameters  is 
virtually  complete  from  the  broad  wavelength  distribution  alone. 


6.  CONCLUSION 

I  have  documented  the  operation  of  a  subroutine, 
PCDIST,  which  may  be  called  from  a  Fortran  program  to  calculate 
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PROGRAM 

VERSION 


EX3.FOR 

FrI  Feb  2618:28:28  1993 


SOURCE  CODE 


PAGE  1 


PROGRAM  EX3 

Main  progron  exoHpIe  3. 


REAL*8  PAR{5000.7).  PF(0:3600).  CSCA,  CABS.  OELTHA,  WAV,  A 
INTEGER  I.  J.  K,  NROUS.  NANG 

OELTHA  -  1 
NANG  -  181 
NROWS  °  3111 

00  10  N1.61 

WAV  -  .20  ♦  REAL(I-1)».03 
00  20  K«1.51 
A  -  1.0  ♦  REAL(K-1)*.04 
PAR((I-l)*51+K.l)  •  3.897*A/WAV 

PAR{(I-1)*51*I£.2)  •  1.52122  ♦  .01095/WAV**2  -  .00000683/WAV**4 

PAR{(I-1)*51+K.3)  -  0.0 

PAR((I-1)*5UK.4)  -  8 

PAR((I-lj*51+K,5)  »  1.2407 

PAR((I-1)*51*K.6)  *  6 

PAR((I-1)*51+K.7)  -  (1/(WAV**5»(£XP(2.4/WAV)-1))) 
*EXP(-.5*((A-2.0)/.3)**2) 


20  CONTINUE 
10  CONTINUE 


CALL  PCDIST(NROWS.  PAR,  OELTHA.  CSCA.  CABS.  PF) 

0PEN(UNIT*11.  FILE-'PHASFUNC.OAT') 

DO  40  J*0. (NANG-1) 

40  WRITE(11.21)  DELTHA*REAL(J).  PF(J) 

21  F0RMAT(F8.2.  £12.4) 

CLOSE (11) 

STOP 

END 


Figure  8.  Fortran  Code  to  Solve  Example  3 
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the  cross  sections  and  phase  function  of  a  cloud  of  irregular 
particles,  using  the  Pollack  and  Cuzzi  approximations.  The 
subroutine  interacts  with  the  programmer's  main  program  through 
structure  called  the  parameter  array,  which  specifies  the 
particles'  properties  through  rules  explained  in  this  report. 
Finally,  to  elucidate  the  use  of  the  parameter  array,  we  have 
given  some  bare-bones  examples  of  main  programs  calling  PCDIST. 
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